--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Commit 9fe7632e8f8285c8c0a17a9966371d0666b05a0c
Parents : ecb5f0c
Author : Mark Qvist <mark@unsigned.io>
Date : 2023-10-25T02:57:28+02:00
Added angle-to-horizon and radio horizon calculations
Changes
Diff
diff --git a/sbapp/sideband/geo.py b/sbapp/sideband/geo.py
index 09326072..c4ca0d53 100644
--- a/sbapp/sideband/geo.py
+++ b/sbapp/sideband/geo.py
@@ -223,42 +223,78 @@ def orthodromic_distance(c1, c2, ellipsoid=True):
else:
return spherical_distance(c1, c2)
-# def tests():
-# import RNS
-# import numpy as np
-# from geographiclib.geodesic import Geodesic
-# geod = Geodesic.WGS84
-# coords = [
-# [(51.2308, 4.38703, 0.0), (47.699437, 9.268651, 0.0)],
-# [(51.230800, 4.38703, 0.0), (51.230801, 4.38703, 0.0)],
-# [(35.3524, 135.0302, 100), (35.3532,135.0305, 500)],
-# [(57.758793, 22.605194, 0.0), (43.048838, -9.241343, 0.0)],
-# [(0.0, 0.0, 0.0), (0.0, 0.0, 0.0)],
-# [(-90.0, 0.0, 0.0), (90.0, 0.0, 0.0)],
-# [(-90.0, 0.0, 0.0), (78.0, 0.0, 0.0)],
-# [(0.0, 0.0, 0.0), (0.5, 179.5, 0.0)],
-# [(0.7, 0.0, 0.0), (0.0, -180.0, 0.0)],
-# ]
-# for cs in coords:
-# c1 = cs[0]; c2 = cs[1]
-# print("Testing: "+str(c1)+" -> "+str(c2))
-# us = time.time()
-# ld = c1+c2; g = geod.Inverse(c1[0], c1[1], c2[0], c2[1])
-# print("Lib computed in "+str(round((time.time()-us)*1e6, 3))+"us")
-# us = time.time()
-# eld = orthodromic_distance(c1,c2,ellipsoid=True)
-# if eld:
-# print("Own computed in "+str(round((time.time()-us)*1e6, 3))+"us")
-# else:
-# print("Own timed out in "+str(round((time.time()-us)*1e6, 3))+"us")
-# ed_own = euclidian_distance(c1,c2,ellipsoid=True)
-# sd_own = orthodromic_distance(c1,c2,ellipsoid=False)
-# aa = azalt(c1,c2,ellipsoid=True)
-# fac = 1
-# if eld: print("LibDiff = "+RNS.prettydistance(g['s12']-eld)+f" {fac*g['s12']-fac*eld}")
-# print("Spherical = "+RNS.prettydistance(sd_own)+f" {fac*sd_own}")
-# # print("EllipLib = "+RNS.prettydistance(g['s12'])+f" {fac*g['s12']}")
-# if eld: print("Ellipsoid = "+RNS.prettydistance(eld)+f" {fac*eld}")
-# print("Euclidian = "+RNS.prettydistance(ed_own)+f" {fac*ed_own}")
-# print("AzAlt = "+f" {aa[0]} / {aa[1]}")
-# print("")
+def distance_to_horizon(c, ellipsoid=False):
+ if ellipsoid:
+ raise NotImplementedError("Distance to horizon on the ellipsoid is not yet implemented")
+ else:
+ # TODO: This is a only barely functional simplification.
+ # Need to calculate the geodesic distance to the horizon
+ # instead.
+ if len(c) >= 3:
+ r = mean_earth_radius
+ h = c[2]
+ return sqrt(pow((h+r),2) - r*r)
+ else:
+ return None
+
+def angle_to_horizon(c, ellipsoid=False):
+ if ellipsoid:
+ raise NotImplementedError("Angle to horizon on the ellipsoid is not yet implemented")
+ else:
+ r = mean_earth_radius
+ h = c[2]
+ return degrees(-acos(r/(r+h)))
+
+def radio_horizon(c1, c2, ellipsoid=False):
+ # dr = 4.12*(√h1 + √h2)
+ if ellipsoid:
+ raise NotImplementedError("Radio horizon on the ellipsoid is not yet implemented")
+ else:
+ h1 = c1[2]
+ h2 = c2[2]
+ ed = euclidian_distance(c1,c2)
+ rh1 = 1e3*4.12*(sqrt(h1))
+ rh2 = 1e3*4.12*(sqrt(h2))
+ rhc = 1e3*4.12*(sqrt(h1) + sqrt(h2))
+ return (rh1, rh2, rhc, rhc > ed)
+
+def tests():
+ import RNS
+ import numpy as np
+ from geographiclib.geodesic import Geodesic
+ geod = Geodesic.WGS84
+ coords = [
+ [(51.2308, 4.38703, 0.0), (47.699437, 9.268651, 0.0)],
+ [(51.2308, 4.38703, 0.0), (47.699437, 9.268651, 30.0*1e3)],
+ # [(51.230800, 4.38703, 0.0), (51.230801, 4.38703, 0.0)],
+ # [(35.3524, 135.0302, 100), (35.3532,135.0305, 500)],
+ # [(57.758793, 22.605194, 0.0), (43.048838, -9.241343, 0.0)],
+ # [(0.0, 0.0, 0.0), (0.0, 0.0, 0.0)],
+ # [(-90.0, 0.0, 0.0), (90.0, 0.0, 0.0)],
+ # [(-90.0, 0.0, 0.0), (78.0, 0.0, 0.0)],
+ # [(0.0, 0.0, 0.0), (0.5, 179.5, 0.0)],
+ # [(0.7, 0.0, 0.0), (0.0, -180.0, 0.0)],
+ ]
+ for cs in coords:
+ c1 = cs[0]; c2 = cs[1]
+ print("Testing: "+str(c1)+" -> "+str(c2))
+ us = time.time()
+ ld = c1+c2; g = geod.Inverse(c1[0], c1[1], c2[0], c2[1])
+ print("Lib computed in "+str(round((time.time()-us)*1e6, 3))+"us")
+ us = time.time()
+ eld = orthodromic_distance(c1,c2,ellipsoid=True)
+ if eld:
+ print("Own computed in "+str(round((time.time()-us)*1e6, 3))+"us")
+ else:
+ print("Own timed out in "+str(round((time.time()-us)*1e6, 3))+"us")
+ ed_own = euclidian_distance(c1,c2,ellipsoid=True)
+ sd_own = orthodromic_distance(c1,c2,ellipsoid=False)
+ aa = azalt(c1,c2,ellipsoid=True)
+ fac = 1
+ if eld: print("LibDiff = "+RNS.prettydistance(g['s12']-eld)+f" {fac*g['s12']-fac*eld}")
+ print("Spherical = "+RNS.prettydistance(sd_own)+f" {fac*sd_own}")
+ # print("EllipLib = "+RNS.prettydistance(g['s12'])+f" {fac*g['s12']}")
+ if eld: print("Ellipsoid = "+RNS.prettydistance(eld)+f" {fac*eld}")
+ print("Euclidian = "+RNS.prettydistance(ed_own)+f" {fac*ed_own}")
+ print("AzAlt = "+f" {aa[0]} / {aa[1]}")
+ print("")
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────